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The exact multiple sphere superposition method is used to calculate the coherent and 
incoherent contributions to the ensemble-averaged electric field amplitude and Poynting 
vector in systems of randomly positioned nonabsorbing spherical particles. The target 
systems consist of cylindrical volumes, with radius several times larger than length, 
containing spheres with positional configurations generated by a Monte Carlo sampling 
method. Spatially dependent values for coherent electric field amplitude, coherent 
energy flux, and diffuse energy flux, are calculated by averaging of exact local field and 
flux values over multiple configurations and over spatially independent directions for 
fixed target geometry, sphere properties, and sphere volume fraction. Our results reveal 
exponential attenuation of the coherent field and the coherent energy flux inside the 
particulate layer and thereby further corroborate the general methodology of the 
microphysical radiative transfer theory. An effective medium model based on plane 
wave transmission and reflection by a plane layer is used to model the dependence of the 
coherent electric field on particle packing density. The effective attenuation coefficient of 
the random medium, computed from the direct simulations, is found to agree closely 
with effective medium theories and with measurements. In addition, the simulation 
results reveal the presence of a counter-propagating component to the coherent field, 
which arises due to the internal reflection of the main coherent field component by the 
target boundary. The characteristics of the diffuse flux are compared to, and found to be 
consistent with, a model based on the diffusion approximation of the radiative transfer 
theory. 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

Modern, parallel-based computational hardware, 
coupled with superposition methods for calculating the 
net scattered field by a collection of particles, have made 
feasible the use of direct simulation strategies to investi- 
gate the characteristics of electromagnetic energy trans- 
port in dense particulate systems. A direct simulation, 
by definition, provides an exact description of the 
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electromagnetic field both within and external to a target 
containing a large yet finite number of particles, and 
which is excited by an external source of radiation. In 
this sense, results generated by direct simulations can be 
considered a benchmark for evaluation of analytical or 
phenomenological theories to describe the propagation of 
electromagnetic energy in discretely inhomogeneous 
media. As a case in point, we have recently used the 
multiple sphere T matrix (MSTM) code to calculate the 
scattering matrices of targets containing up to several 
thousand spheres, with the individual spheres having size 
parameters up to 4 [1-3]. The targets, in these calcula- 
tions, consisted of spherical volumes, with the spheres 
randomly distributed within the volume with a set overall 
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volume fraction; volume fractions used in the calculations 
ranged from 0.01 to 0.4. Since the number of spheres in 
the targets was large, we assumed that the configuration- 
averaged scattering matrix (i.e., that obtained from aver- 
aging scattering matrices, each obtained for a fixed 
realization of the sphere properties and a fixed incident 
direction) could be approximated by the orientation 
averaged scattering matrix for a single realization; this 
allowed us to exploit the analytical orientation averaging 
properties of the target T matrix as calculated by MSTM. 
Results from these calculations have provided definitive 
evidence to support microphysical theories of radiative 
transfer and coherent backscattering [4], 

The objective of this work is to present a similar 
comparison between numerically exact computer simula- 
tions and analytical theory. The phenomena under exam- 
ination, in this case, are the attenuation of the coherent 
field, and the propagation of the coherent and diffuse 
components to the time-averaged Poynting vector in a 
particulate medium, due to scattering by the particles. 
The simulation procedure is relatively straightforward: 
local electric and magnetic field amplitudes will be 
calculated for a statistically representative sequence of 
randomly generated targets of spheres, with each target 
configuration corresponding to the same sphere proper- 
ties and average concentration. The local coherent field 
and time-averaged Poynting vector are then obtained by 
assuming ergodicity and averaging over the sequence of 
configurations. An important difference in the procedures 
used here, as opposed to our previous investigations [2,3], 
is that they will not employ analytical T matrix averaging: 
our focus is on the local field amplitudes within the 
medium as opposed to the far-held scattering behavior. 

A well established theoretical framework has been 
developed to describe coherent held attenuation in dis- 
cretely inhomogeneous media. On the most basic level is 
the standard radiative transfer theory, which represents 
the dilute concentration limit [2,4]. The imaginary part of 
the effective propagation constant in the medium, for this 
approximation, is obtained from the product of the 
particle number density and average extinction cross 
section. 

The radiative transfer theory is based on the far-held 
version of the Foldy-Lax equations [4] and will fail for 
wavelength-sized particles with concentrations appro- 
aching and exceeding 0.1. For such cases comprehensive 
effective medium (EM) formulations have been devel- 
oped, most notably by Ishimaru and the Varadans [5-7], 
Such theories begin with the same basic superposition 
model employed in the direct simulation, for which the 
exciting Held at a particular sphere is given by the 
incident (i.e., externally exciting) held and the sum of 
helds scattered from every other sphere in the system. 
Unlike a direct simulation - which solves the superposi- 
tion equations for each sphere in the system and then 
averages over multiple configurations - EM formulations 
attempt to average the superposition equations analyti- 
cally, with subsequent derivation of relations for the 
effective propagation constant. Performing this averaging 
requires knowledge of the pair correlation function of the 
particles in the random medium, and it also invokes 


simplifying assumptions regarding the high-order corre- 
lations among the particle positions and scattered Helds, 
e.g., the quasi-crystalline approximation (QCA) [8-10]. 

Both computational and experimental evidence has 
been collected to validate EM theories. Tsang et al. [11] 
inferred the attenuation rate from direct simulations of 
far-held scattering by a particulate volume and obtained 
good correspondence with theory, although the simula- 
tions were limited to spheres with a relatively small size 
parameter of 0.2. Analytical estimates of EM theories have 
been shown to accurately predict laboratory measure- 
ments of extinction in random suspensions of monodis- 
perse spheres having size parameters of order unity or 
greater [12]. 

The investigations in [11,12] were focussed entirely on 
the effective propagation constant (or, equivalently, effec- 
tive refractive index) of the spherical-particle system. 
As predicted by theory, this quantity depends solely on 
the sphere properties, concentration, and pair correlation 
function; it is not a function of the geometry of the system 
containing the particles. Unlike Ref. [11], the objective of 
this work is to examine how the target geometry affects 
the coherent Held within the target. In particular, we will 
employ targets which model a Hxed-thickness slab of 
particles. In an EM model, such a system would represent 
a plane layer. It is well known that interference of the 
forward and backward-propagating plane waves in a 
homogeneous plane layer can strongly affect the reflec- 
tance and transmission of the layer. We will show that the 
coherent Held, in a slab of random spherical particles, can 
also cause interference effects. 

An additional objective of this work is to examine the 
character of electromagnetic energy transport in the 
random medium. In particular, the direct simulations, 
combined with configurational averaging, will be used 
to calculate the coherent and diffusive components to the 
time-averaged Poynting vector in the medium. We will 
apply the diffusion approximation of the radiative trans- 
fer theory to model the diffuse flux. 

2. Methodology 

2.2. Basic definitions 

We assume that the time dependence of the electric 
and magnetic helds is harmonic and described, in the 
complex-held representation, by the simple complex 
exponential exp(-icot), where co is the angular frequency, 
t is time, and i = v^T. In other words, we assume that the 
complex electric and magnetic Helds can be factorized as 
E(r,t) = exp(-icot)E(r) and H(r,t) = exp(-icot)H(r) respec- 
tively, where r is the position (radius) vector, while the 
actual real-valued helds are obtained by talking the real 
part of the respective complex Helds. The complex held 
amplitudes E(r) and H(r) are constant if the scattering 
target is Hxed, but vary with time implicitly if the 
scattering object undergoes temporal changes. In what 
follows, we will assume that temporal huctuations of E(r) 
and H(r) caused by changes in the scattering object occur 
much more slowly than the harmonic oscillations of the 
factor exp(-icot). 


D.W. Mackowski, M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 123 (2013) 103-112 


105 


It is well known (see, e.g., Ref. [4]) that the instanta- 
neous direction and rate of the electromagnetic energy 
flow at r are described by the Poynting vector 
S(r,t) = Re[E(r,f)] x Re[H(r,t)]. Averaging the Poynting vec- 
tor over a time interval capturing a sufficiently large 
number of time-harmonic oscillations of the electromag- 
netic field yields <S(r,t)> =Re[S(r)], where the complex 
Poynting vector is defined by S(r) = (l/2)E(r) x H*(r) and 
the asterisk denotes a complex-conjugate value. 

Let us consider the scattering object in the form of a 
volume of discrete random medium. The statistical ran- 
domness of this object is caused by chaotic movements of 
the participating particles resulting in relatively slow 
temporal fluctuations of the complex field amplitudes 
E(r) and H(r). At each given moment these amplitudes can 
be computed by solving the frequency-domain Maxwell 
equations for the corresponding instantaneous multi- 
particle configuration. In the majority of practical applica- 
tions, one needs to calculate the average of a specific 
optical observable such as S(r,t) over a time interval 
typically required to take an actual measurement. In what 
follows, we will assume that the scattering object is 
ergodic (e.g., Ref. [4]) and will replace averaging over 
time by averaging over a large statistically representative 
set of randomly created realizations of the multi-particle 
object. In other words, we will assume that 
C S( r, t) » = Re[«S(r)>] = S(r) where the overbar denotes 
ensemble averaging and C . . . > denotes averaging over a 
time interval much longer than the typical period of 
random fluctuations caused by chaotic movements of 
the particles. 

For the sake of brevity, we will often refer to E(r) and 
H(r) as the electric and magnetic field, respectively, which 
is a common practice. It should be remembered, however, 
that these quantities represent only the time-independent 
(or slowly fluctuating) amplitudes of the actual time- 
harmonic physical fields. 

2.2. Target generation 

The targets used in the simulations consisted of a 
cylindrical volume of radius R and axial length L. Within 
this volume are placed the origins of N s nonoverlapping 
spheres, each with identical radius a and refractive index 
m. The axis of the cylinder is taken to be the z axis in a 
global coordinate system. The quantities N s and L were 
varied in the simulations, whereas R was set to a value 
sufficiently large, relative to both the sphere radius a and 
L, so that the average Poynting vector, calculated in 
regions close to the target axis, was unaffected by the 
lateral boundary of the target. A maximum value of 
L=10a was used in the simulations, and for this length 
a target radius of R = 30 a to 40a was sufficient to simulate 
plane-layer conditions. 

For a set R and L, the positions of the N s spheres in 
the target were generated via a Monte Carlo sampling 
procedure. Specifically, a trial sample of the cylindrical 
coordinates p jt </>,-, z,- of sphere i in the set is obtained 
according to 

Pi = RTZ 1 /2 , <pi = 2 tcTZ, Zj = Zj.i + TZ Az, ( 1 ) 


with 



( 2 ) 


in which TZ denotes a random number uniformly distrib- 
uted between 0 and 1 (each occurrence of TZ is implied to 
be a new sample), and / is the average volume fraction of 
the spheres. This procedure automatically constrains the 
volume fraction of the spheres in the target volume to be/. 
The trial obtained from Eq. (1) is now checked against the 
positions for sphere j — 1 , i—2, . . ., and the trial is accepted 
if no overlap occurs between spheres and is rejected and 
resampled if overlap does occur. A starting value of Zj = 0 
is arbitrarily set, and the procedure is carried out to some 
value N samt> ?z 1.5-2 times N s , where N s denotes the 
desired number of spheres in the target. The actual 
positions of the target spheres are then chosen as those 
of the topmost N s samples in the set; the rationale here is 
that any pair-correlation artifacts in sphere positions 
occurring at startup will be washed out in the topmost 
samples. 

Fig. 1 shows a typical target configuration used in the 
simulations. The target has dimensions R = 40a, L=10a, 
and contains N s = 6000 spheres; the corresponding 
volume fraction is /= 0.5. 


2.3. Calculation procedure 

The Multiple Sphere T Matrix (MSTM) code [1,13] was 
used to calculate the electric E(r), magnetic H(r), and 
Poynting S(r) vector fields in and around a fixed multi- 
particle target for a sequence of sampled targets. Para- 
meters held constant in the target sequence were the 
sphere size parameter and refractive index, the target 
radius R and thickness L, and the number of spheres N s ; 
the latter constraint is equivalent to fixing the average 
volume fraction of the spheres in the target. The target 
was excited with a plane wave, propagating in the x-z 
plane at an angle jf relative to the positive direction of the 
z axis and linearly polarized either parallel or perpendi- 
cular to the x-z plane. 

The formulation employs a superposition representa- 
tion of the electric field exterior to the spheres, that being 

Ns 

Eew(r) = E inc (r) + E sra (r) = E inc (r) + J2 E S ca,i(r) (3) 

1 = 1 



Fig. 1 . Random sphere target. 
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The scattered field, from a particular sphere, is described 
by outgoing vector spherical wave functions (VSWFs) 
centered about the origin of the sphere, and the incident 
field is described by a regular VSWF expansion centered 
about a common target origin: 

Ei„c(r)= E E/n.np N ^p(k 0 r) (4) 

n = im = -np = i 



(5) 


In the above, k 0 is the wavenumber in the external 
medium, r f is the position vector of the center of the ith 
sphere, and N mnp denotes the VSWF of either type 1 
(regular) or 3 (outgoing), of order n, degree m, and mode 
p= 1 (TM) or 2 (TE). The truncation orders for the spheres, 
L s , is set by inspection, i.e„ by comparison of results 
computed with successively higher L s . For the conditions 
examined here, L s was not dependent on the sphere 
concentration, and the required number of orders would 
depend solely on the size parameter and refractive index 
of the spheres. The truncation order for the incident field, 
L t , will typically scale with the dimensionless distance 
k 0 1 r | of the evaluation point from the target origin, and 
the incident field coefficients f mnp will depend on the 
direction and propagation state of the incident field. 

The translation theorem for VSWFs, combined with the 
Mie relations for the sphere, results in a system of 
equations for the scattering coefficients a‘ mnp : 


1 

5np 


a 


mnp 


N s L s l 2 

V s H l ~j d —f P ik o(x« sin p+Zi cos p) 

/ , / v / v / , n mnpklq u klq ~Jmnp c 


i = } l=\k = —lq=\ 


( 6 ) 


where a np are the Mie coefficients, which are functions of 
the sphere size parameter k a and refractive index m, and 
H'~i is the outgoing translation matrix, which depends 
solely on the distance and direction of translation from 
origins j to i. 

Eq. (6) represents 2N s L s (L s +2) complex valued equa- 
tions for the set of scattering coefficients. Iterative meth- 
ods, combined with techniques to accelerate the 
translation operations, are employed in the MSTM code 
to generate a solution. All relevant microscopic and 
macroscopic electromagnetic properties of the target, 
i.e., external and internal field amplitudes, cross sections, 
scattering matrix elements, can be obtained analytically 
from the solution, and the predicted values are exact to 
the truncation error of the VSWF expansions. Of particular 
relevance here is the calculation of the complex electric 
and magnetic field vectors at arbitrary points. This is 
accomplished using the superposition of scattered and 
incident fields in Eqs. (3)-(5) for points exterior to the 
spheres, while at interior points, say within sphere i, the 
electric field is given by an expansion of regular VSWFs, 
evaluated in the sphere medium 


L s n 2 

E inU (r)=J2 J2 H c mnp N mnp(mko(r-r,)) 
n = lm = -np=l 


in which the internal field coefficients cj„ np are obtained 
from the scattering coefficients aj,, np via the Mie relations. 


2.4. Configurational averaging 

The basic strategy of the calculations was to calculate 
and record electric and magnetic field amplitudes on a set 
of fixed sample points r p , p = l,2, ..., for a sequence of 
randomly generated sphere configurations, with each 
configuration corresponding to a set target geometry, 
sphere properties, and sphere volume fraction. The sam- 
ple points are placed on a grid within a rectangular 
sample volume centered on the target axis. This sample 
volume, which is depicted in Fig. 2, has dimensions of 
W xWx H, where the lateral length W is typically around 
1/4 to 1/2 of the target radius R, and the axial height H is 
set so the volume extends both above and below the top 
and bottom surfaces of the target. 

The so-called coherent electric field, at a fixed sample 
point r p , is defined as the average of the electric field 
amplitude E(r) over a sufficiently long period of time. 
Assuming ergodicity, it can be obtained by averaging the 
calculated field amplitude values at the point over the 
sequence of random configurations, i.e., 



where c k denotes the parameters for configuration k and 
N c is the number of sampled configurations. As explained 
in Ref. [14], the coherent electric field is not a real 
physical field since it is obtained by averaging the com- 
plex amplitude of the electric field rather than the electric 
field itself. The practical usefulness of this purely math- 
ematical vector field and its magnetic analog will be clear 
from the following discussion. 

Since the incident field propagates in the x-z plane - 
and, for a target of infinite lateral extent (R-»oo), the 
coherent field would be independent of y - we would 
further refine the coherent field estimate by averaging 
over the sample points in the y direction for a fixed x p , z p , 



(7) 


Fig. 2. Target and sampling configuration. 
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via 


E(^p.Zp) 


1 

Np 


E 


E (x p ,y p ,Zp) 


(9) 


where N p denotes the number of points in the y direction. 
For the specific case of normal incidence (/? = 0), the 
coherent field would depend solely on z, and averaging 
was performed over both the x and y directed sample 
points. 

To avoid complicating the subsequent analysis, we will 
not distinguish between the normal and oblique inci- 
dence cases and their associated strategies for coherent 
field calculation: the coherent field at a point r will be 
denoted simply as E(r), and it is implicitly understood 
that calculation of this quantity involves both configura- 
tional and spatial averaging. Given this, the electric field 
amplitude at a specific sample point r p , corresponding to 
a specific configuration c, can now be decomposed into 
the coherent and incoherent parts [5,6], via 

E(r p ,c) = E(r p ) +E'(r p ,c) (10) 

where the ensemble average of E'(r p ,c) is zero. An analo- 
gous decomposition holds for the magnetic field at point 
r p and configuration c. 

The short-time average of the real-valued Poynting 
vector, at a sample point r p and for a specific target 
configuration c, is given by 

< S(r p , c,t)> =iRe[E(r p ,c) x H*(r p ,c)] (11) 

The long-time average of this quantity at a point r p can be 
calculated per the strategy used to calculate the coherent 
field, i.e., assuming ergodicity and averaging over config- 
urations and over the spatially independent direction in 
the sample volume. By the use of Eq. (10), the Poynting 
vector, averaged over a sufficiently long period of time, 
can be decomposed into a coherent part and a diffuse part 
[5,6]: 

< S( r p , c, t) > = S( r p ) = S coh (r p ) + S d ( r p ) (12) 

in which 

Scoh(r P ) = 2 Re[E(r p ) x H (r p )] (13) 


Sd(r P ) = \ Re[E'(r p ,c) x H'*(r p ,c)] (14) 

The quantities that can be calculated directly, on the 
set of grid points r p , are the coherent fields E(r p ) and 
H(r p ), the average Poynting vector S(r p ), and the coherent 
flux S coh (r p ). The scattering flux is obtained from Eq. (12). 


3. Results 

3.1. General trends 

We begin by examining the general characteristics of 
the coherent field and the coherent and diffuse fluxes for 
targets exposed to oblique-incidence plane waves. The 
target, for this case, has a dimensionless radius k 0 R = 60, 
dimensionless thickness k 0 I = 15, and contains N s = 2250 
spheres, each with size parameter k 0 a = 2 and refractive 
index m = 1.48 (characteristic of quartz in the visible 
wavelengths); the average volume fraction for this case 
is /= 0.25. Three orders were used to represent the 
scattered fields from the spheres, and the interaction 
equations in Eq. (6) contained 67,500 unknowns for each 
incident state and target configuration. The maximum 
order required to evaluate the incident field at points 
within the sample volume, via Eq. (4), was approximately 
L t 60. The MSTM code was run on the Auburn Univer- 
sity College of Engineering compute cluster, using 64 
processors, and complete solution for a set target config- 
uration and incident angle (corresponding to two solu- 
tions to Eq. (6) for two mutually perpendicular incident 
polarizations) required around 15 min. A total of 60 
random configurations were used to compute the coher- 
ent fields and extinction and scattering fluxes. 

Fig. 3 contrasts the structure of the local electric field 
amplitude for a particular configuration vs. that for the 
coherent field following configurational and spatial aver- 
aging. Both plots show levels of Re y ■ E in the x-z plane. The 
incident field has /? = 30°, and is polarized perpendicular to 
the plane. Axes coordinates are k 0 x and k 0 z, with the center 
of the target corresponding to the origin, and the target is 
contained within k 0 z = + 7.5. The coherent field on the right 
has been averaged over 60 configurations and over 21 planes 
in the y direction evenly spaced between k 0 y = ± 20. 



Fig. 3. Local value of Re(y • E) (left), and coherent field Re(y • E). 
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Fig. 4. Vector field of S in x-z plane (top left), averaged (top right), and coherent (bottom left) and diffuse (bottom right) components of S. Shading 
represents magnitude of vector. 


The coherent field in Fig. 3 contains characteristics 
associated with the reflection and transmission of a plane 
wave by a homogeneous layer. In particular, a refraction 
and retardation of the wave is evident as the field 
propagates through the target medium. This behavior is 
more apparent when the average Poynting vector is 
decomposed into coherent and diffuse parts. Shown in 
Fig. 4 are stream plots of the local Poynting vector in the 
y = 0, x-z plane for a fixed configuration, along with total 
average Poynting vector, and the coherent and diffuse 
contributions to the total average, for /? = 30° incidence 
and perpendicular polarization. 

The two top plots in Fig. 4 display the directional flow 
of the total energy transfer in and about the target, with 
the left being a “snapshot” for a particular configuration, 
and the right showing the average flow pattern. The 
average flux entering the system is primarily in the 
direction of the incident beam, and it is steered toward 
the normal direction during the transport across the 
target. A simple radiative transport interpretation of this 
effect is that scattering by the particles results in an 
effective retardation of the energy flow, and the effect of 
this is to shift the overall direction of the flow toward the 
normal directions. 

A more concrete interpretation, in our view, makes use 
of the principle of energy conservation. The spheres, for 
this specific example calculation, are nonabsorbing, and 
for such a situation S must be divergence free at all points. 
Furthermore, for a target of infinite lateral extent, this 
condition would imply z S = constant at all axial posi- 
tions. Scattering by the spheres (and absorption, if pre- 
sent) act to extinguish the coherent field - as observed in 
Fig. 3 - and this results in a sink in the coherent flux per 
the definition in Eq. (13). Due to overall conservation of 
energy, this sink must reappear as a source in the diffuse 
flux. The two lower stream plots in Fig. 4 illustrate this 
interaction between coherent field destruction and diffuse 
flux (or, equivalently, incoherent field) formation. 
The coherent flux vector field behaves analogously to a 
plane wave interacting with a absorbing homogeneous 


layer, in which the effective absorption describes the 
destruction of the coherent field by scattering. Likewise, 
the diffuse flux is created within the target, and emerges 
from both the top and bottom surfaces. 


3.2. Effective medium model of coherent field 


We will now focus attention on the extinction of the 
coherent field within the medium, with the objective of 
identifying, via direct simulations coupled with modeling, 
the effective refractive index of the target. To simplify the 
presentation, we will restrict the examination to the 
normal incidence case (fl = 0). 

The effective medium theories of [5,7] were based on a 
half-space model of the inhomogeneous medium - for 
which the regions z < 0 and z > 0 consisted of free space 
and a random distribution of particles, respectively - that 
was excited by an x polarized, unit amplitude plane wave 
propagating in the +z direction. The coherent field, for 
this case, is a plane wave characterized by an effective 
propagation constant, i.e., 

E C oh(z) = xe ik " z (15) 


The most straightforward estimation of the effective 
refractive index m e = k e /k 0 , from the simulation results, 
would therefore be obtained from 


m e 


1 l n /x • E CO h(L)^ 

ik of yx ■ E coh (0) j 


(16) 


where E co /,(0) and E coh (L) are the calculated complex 
coherent field at the target boundaries. 

The half space model of Eq. (15) neglects the existence 
of a downward-propagating component in the coherent 
field, which could arise from a reflection of the upward- 
propagating component from the top boundary. To exam- 
ine such effects, we performed a more detailed modeling 
of the coherent field based on a homogeneous layer 
(or film) of thickness L e and refractive index m e . The model 
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formula for E coh (z ) is given by 


( e ik ° z 


+ re 


-ikoz 


Ecofi(z) =x< ae' Kz + be lk ' z , 
I te ik ° z , 

with 


z< 0 
0 < z < L e 
z>L e 


(17) 


-[(m 3 -l)(l-e 2im " 1 ')] 

(18) 

2(m e + 1) 
D 

(19) 

2(m e -l)e 2im ' L ' 

D 

(20) 

4m f ,e l(m '~ I) re 

D 

(21) 

• 2im ’Li’(-i-rn e ) 2 _(1+m e ) 2 

(22) 


b= - 


t=- 


and 


The coherent field model has three parameters, being 
the real and imaginary parts of m e and the effective 
thickness L e . We keep the thickness L e a parameter - as 
opposed to setting it equal to thickness L used to generate 
the targets - because the target does not have sharp, 
distinct boundaries. That is, since all sphere origins must 
lie between z=0 and L, there will be a region, of thickness 
equal to the sphere radius, over which the medium 
transitions from inhomogeneous to homogeneous. Values 
of m e and L e , for a given set of coherent field simulation 
data, were obtained by a nonlinear least square error fit of 
Eq. (17) to the data. To perform this fit, the point z = L e /2 
in Eq. (17) was made to coincide with the target center in 
the simulation data. 

Simulations were performed using spheres with 
k 0 a = 2.645 and m = 1.194, which correspond to the latex 
spheres-in-water experimental measurements conducted 
in [15] and examined theoretically in [16]. The target 
thickness and radius were k 0 E = 26.45 and k 0 R = 79.35 
(i.e., 10 and 30 sphere radii), and volume fractions ranged 
from f= 0.01 to 0.5. 


Results appear in Fig. 5, which show the coherent field 
phase and magnitude as a function of position, for four 
values of volume fraction /ranging from 0.025 to 0.5. The 
solid lines are the simulation results, and the dashed lines 
show the predictions of Eq. (17) using the fitted m e and 
k 0 i e values for the data set. The simulation results and 
model predictions are presented so that z=0 represents 
the lower boundary of the slab, and the vertical colored 
lines, in the neighborhood of k 0 z= 25-30, mark the fitted 
value of k 0 L e for the particular set. The phase results 
indicate the increasing retardation of the coherent wave 
as concentration increases; results for/=0.025 are not 
shown in this plot because they are negligibly different 
than those for /= 0.079. The magnitude results, which are 
plotted on a log scale, show an exponential decrease in 
amplitude. Of particular interest are the oscillations in 
E(z)j seen for the relatively large / values. This can be 
interpreted as an effect of interference between the 
downward and upwards propagating wave components 
that are present in a finite-thickness layer. The model 
results overestimate the amplitude of the downward 
component - likely because of the representation of the 
boundaries as sharp interfaces - yet the simulation results 
do reveal that a downward-propagating component is 
present in the coherent field. 

The specific values of sphere size parameter and 
refractive index, used in the results of Fig. 5, were chosen 
to allow a comparison of the predicted attenuation 
coefficient a t = 2k 0 lm(m e ) of the medium with that 
derived from effective medium theory, and measured by 
experiment, as described in [16]. A useful indicator of the 
packing density effects on attenuation is the ratio of the 
exact result to that predicted from radiative transfer 
theory, given by 

ct t 4 Im(m e )k 0 a 

1 3/C rat /(47ta 3 ) 3 fQ ext 1 

where C ex[ and CW are the extinction cross section and 
extinction efficiency of the sphere. For sufficiently small / 
this ratio should go to unity. Results are given in Fig. 6, 
which show attenuation ratios as a function of volume 
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Fig. 5. Re x E CO h(z) (left) and \E coh (z)\ (right) vs. k 0 z, for k 0 a = 2.645 and m = 1.194 and for different volume fractions. Solid line: simulations, dashed 
lines: model fit to Eq. (17). Colored vertical lines denote fit values of kL e . 
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Volume fraction/ 


Fig. 6. Attenuation ratio y (Eq. (23)) vs. volume fraction/, for simulation 
results predicted by the layer model (Eq. (17)) and the half-space model 
(Eq. (16)), along with experimental and theoretical results of [16]. 
Sphere and target properties are the same as Fig. 4. 

fraction /. Two sets of simulation results are presented, 
corresponding to oc t derived the simple plane wave model 
of Eq. ( 1 6) and the layer model of Eq. ( 1 7 ). Also appearing 
are the theoretical and experimental results reported in 
[16]. 

The results in Fig. 6 show that, for the specific condi- 
tions examined here, the attenuation coefficient can be 
accurately predicted by treating the target as a half 
space — in that the values of a t derived from Eq. (16) are 
not significantly different than those obtained from a 
more comprehensive model. Furthermore, the simulation 
results are highly consistent with both effective medium 
theory and experiment, especially for / greater than 
around 0.1. For smaller volume fractions the simulation 
results show a small degree of scatter around the asymp- 
totic limit of unity. We suspect that this is an artifact of 
the averaging process, as a smaller number of spheres in 
the target - corresponding to smaller / - leads to a 
relatively larger standard deviation in the configuration- 
dependent properties of the target. Consequently, a num- 
ber of configuration samples greater than the 100 used 
here for each point may be required to accurately resolve 
the coherent field properties, yet, at the same time, the 
simple, independent scattering model works quite well in 
this limit. 


3.3. The diffuse flux 

Per the discussion in Section 3.1, the divergence of the 
average Poynting vector S(r) is zero for the nonabsorbing 
spheres examined here. The destruction of the coherent 
flux S co /,(r) by “effective" absorption - as manifested in a 
nonzero imaginary part of the effective refractive index 
m e - must be balanced by the creation of the diffusive flux 
SdW- 


Ishimaru et al. [17] applied the diffusion approxima- 
tion to the radiative transport equation to model the 
diffuse intensity and flux in systems containing relatively 
large densities of scattering particles. It should be empha- 
sized that the intensity - being, per the common radiative 
transport definition, a scalar quantity describing the flow 
of energy in a particular direction and contained within a 
differential solid angle - cannot be extracted from the 
simulations results [18]. That is, the only physically 
meaningful quantities obtained from the simulations are 
the electric and magnetic field amplitudes and the Poynt- 
ing vector. In this sense, a comparison of a radiative 
transport equation solution for the diffuse intensity with 
a direct simulation, based on Maxwell’s equations, is 
impossible. However, a radiative transport prediction of 
the diffuse flux, which is obtained from an integral of the 
diffuse intensity, can be compared with the simulation 
results. 

To do this, we adapt the diffusion approximation 
developed in Ref. [17] to the conditions of the simula- 
tions. In particular, under normal incidence conditions, 
the diffuse flux will be parallel to the z axis and will be a 
function solely of z. For the nonabsorbing conditions 
examined here (albedo of unity), the diffuse flux will be 
related to the average diffuse intensity U(z ) via, 


Sd,z(Z) = 


1-gi 



47t dU{z)\ 
3 x t dz ) 


(24) 


where S d z (z) = z ■ S d (z), f 0 is the intensity of the incident 
field ( = 1 /2 for the unit-amplitude incident electric field), 
and gi is the dimensionless scattering asymmetry para- 
meter of the particles (i.e., a characteristic of the fraction 
of energy single-scattered by the particles into the for- 
ward directions). The average diffuse intensity, for non- 
absorbing conditions, satisfies 

dhm = _3i^ e _ XtZ 

dz 2 4 7t 


with the boundary conditions, 

U(0) = -^-S^(0) (26) 


U(L) — (27) 

where the fluxes at the boundaries, S dz (0) and S dz (L), are 
related to U via Eq. (24). The solution for the diffuse flux is 
obtained as 

/ 5-e~ x < L \ 

1281 

To compare the model with simulation results, we set 
a t = 2k 0 Im(m e ) and L = L e , where m e and L e correspond to 
the values obtained, for a specific volume fraction f by 
fitting the coherent field to the plane layer model. 
We then determined a value of the asymmetry parameter 
gi, for each f by minimizing the least square error 
between the simulation and diffusion model predictions 
of the diffuse flux. 

Shown in Fig. 7 are the simulation and model predic- 
tions of the diffuse flux, as a function of dimensionless 
position k 0 z. For / up to around 0.2, the diffusion model 
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Fig. 7. Axial component of diffuse flux S dz (z ) vs. dimensionless axial 
position k 0 z, with volume fraction / as a parameter. Solid lines are 
simulation results, dotted are model predictions from Eq. (28). 
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Fig. 8. Asymmetry parameter g lt calculated from fitting of Eq. (28) to 
simulation results, vs. volume fraction /. The solid line is the single 
sphere Mie result. 
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does an excellent job of capturing the magnitude and 
distribution of the diffuse flux as calculated by the 
simulations. The simulation results for /= 0.5, and to a 
lesser extent for/=0.22, show characteristics near z= 0 
that could not be represented by the exponential model of 
Eq. (28). It is interesting to note that the oscillations seen 
in the /= 0.5 data have a wavelength roughly equal to 
twice the sphere size parameter, and this may point to an 
effect due to an emerging lattice structure among the 
spheres. 


The fitted values of asymmetry parameter g, are 
shown in Fig. 8 as a function of /. The values of g ^ are 
fairly constant for /up to around 0.2, and are roughly 0.8 
of the single-sphere Mie value. Undue emphasis should 
not be placed on the correspondence (or lack thereof) of 
the fitted gi results to the Mie value. The largest optical 
depth ot t L e among the simulation results is 0.42 for 
/=0.33, and conditions for/ <0.1 could be justifiably 
described as optically thin and would contradict the 
optically thick assumptions inherent in the diffusion 
approximation. What is more relevant, in regard to the 
results in Fig. 8, is the transition in gj from a constant, at 
low / to a decreasing function of/. This behavior indicates 
an increasing effect, for large/, of pair-correlated positions 
among the particles, and is consistent with the results of 
Ref. [19] obtained using the structure-factor approach 
coupled with the Percus-Yevick approximation [20]. 


4. Summary 

We have demonstrated that numerically exact solutions 
of Maxwell’s equations in multiple sphere systems, com- 
bined with configurational and spatial averaging, can pro- 
vide direct simulations of the coherent field and the 
coherent and diffuse energy fluxes that are entirely consis- 
tent with theoretical formulations of radiative transport in 
random, discretely inhomogeneous media. We should 
emphasize that the performed simulations are quite tract- 
able on modern parallel computer platforms; the total 
computational time involved in this investigation was on 
the order of 10’s of hours. Indeed, we will submit that direct 
simulations have a future role not as simply verifying 
radiative transfer theories, but rather as a practical means 
of calculating radiative transfer in high-concentration parti- 
cle systems such as regoliths, pigment layers, deposits, etc. 
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